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We study the two-dimensional kinetic Ising model below its equilibrium critical temperature, 
subject to a square-wave oscillating external field. We focus on the multi-droplet regime where the 
metastable phase decays through nucleation and growth of many droplets of the stable phase. At 
a critical frequency, the system undergoes a genuine non-equilibrium phase transition, in which the 
symmetry-broken phase corresponds to an asymmetric stationary limit cycle for the time-dependent 
magnetization. We investigate the universal aspects of this dynamic phase transition at various tem- 
peratures and field amplitudes via large-scale Monte Carlo simulations, employing finite-size scaling 
techniques adopted from equilibrium critical phenomena. The critical exponents, the fixed-point 
value of the fourth-order cumulant, and the critical order-parameter distribution all are consis- 
tent with the universality class of the two-dimensional equilibrium Ising model. We also study the 
cross-over from the multi-droplet to the strong-field regime, where the transition disappears. 

PACS numbers; 64.60.Ht, 75.10.Hk, 64.60. Qb, 05.40.-a 



I. INTRODUCTION 



Metastability and hysteresis are widespread phenomena in nature. Ferromagnets are common systems that exhibit 
these behaviors , but there are also numerous other examples ranging from ferroelectrics to electrochemical 
adsorbate layers to liquid crystals A simple model for many of these real systems is the kinetic Ising model 
(in either the spin or the lattice-gas representation). For example, it has been shown to be appropriate for describing 
magnetization dynamics in highly anisotropic single-domain nanoparticles and uniaxial thin films [pd|-p^ . 



(a) 





FIG. 1. Metastable decay in the multi-droplet regime at T=0.8rc for an L=128 square-lattice kinetic Ising system evolving 
under Glauber dynamics. The system is initialized with all spins Si = l, and an applied field, H=—0.3J, is set at t=0. Snapshots 
of the spin configurations are given at (a) t=30 Monte Carlo steps per spin (MCSS), (b) t=60 MCSS, (c) t=74 MCSS. The 
metastable lifetime (the average first-passage time to zero magnetization) at this temperature and field is (r)=74.5 MCSS. 
Stable s=— 1 (metastable s=-|-l) spins are represented by black (white). 
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The system response to a single reversal of the "external field" is fairly well understood |15[]. In sufficiently large 
systems below the equilibrium critical temperature, Tc, the order parameter changes its value through the nucleation 
and growth of many droplets, inside which it has the equilibrium value consistent with the value of the applied field, 
as shown in Fig. ^ This is the multi-droplet regime of phase transformation |^,^. The well-known Avrami's law 
describes this process of homogeneous nucleation followed by growth quite accurately up to the time when the 
growing droplets coalesce and the stable phase becomes the majority phase The intrinsic time scale of the system 
is given by the metastable lifetime, (r), which is defined as the average first-passage time to zero magnetization. It is 
a measure of the time it takes for the system to escape from the metastable region of the free-energy landscape. In 
this paper we will use the magnetic language in which the order parameter is the magnetization, to, and its conjugate 
field is the external magnetic field, H . Analogous interpretations, e.g., using the terms polarization and electric field 
for ferroelectric systems 1^,0, and coverage and chemical potential for adsorption problems [^|j9j, are straightforward. 

It is natural next to ask, "what is the response to an oscillating external field?" The hysteretic behavior in 
ferromagnets has attracted significant experimental interest, mainly focused on the characteristic behavior of the 
hysteresis loop and its area. Its dependence on the field amplitude and frequency has been intensively studied 
and its scaling behavior (power law versus logarithmic) is still under investigation, both experimentally [|ll| )l3| and 
theoretically p^pTf . For the kinetic Ising ferromagnet in two dimensions it has been recently shown ]25|-|27| that 
the true behavior is in fact a crossover, approaching a logarithmic frequency dependence only for extremely low 
frequencies. 

An important aspect of hysteresis in bistable systems, which is the focus of the present paper, is the dynamic 
competition between the two time scales in the system: the half-period of the external field ti/2 (proportional to 
the inverse of the driving frequency) and the metastable lifetime (t). For low frequencies, a complete decay of the 
metastable phase almost always occurs in each half-period, just like it does after a single field reversal. Consequently, 
the time-dependent magnetization reaches a limit cycle which is symmetric about zero [Fig. |^(a)]. For high frequencies, 
however, the system does not have enough time to switch during one half-period, and the symmetry of the hysteresis 
loop is broken. The magnetization then reaches an asymmetric limit cycle [Fig. ||(b)]. Avrami's law (l^Jl^ is a good 
approximation when the majority of the droplets do not overlap. Thus, it can be employed to estimate the time- 
dependent magnetization and the dynamic order parameter (period-averaged magnetization) in the low-frequency 
(see the Appendix) and in the high-frequency [|2^] limits. However, it cannot describe the "critical regime" where ti/2 
becomes comparable to (r) , and which is dominated by coalescing droplets. 
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FIG. 2. Monte Carlo magnetization time series (solid lines) in the presence of a square-wave external field (dashed lines) for 
an L=128 system at T=0.8rc and field amplitude Hq—O.SJ. {Tc is the two-dimensional equilibrium Ising critical temperature.) 
The metastable lifetime at this temperature and field is (r)=74.5 MCSS. (a) Dynamically disordered phase at dimensionless 
half-period 0=ti/2/(''")= 2.7. (b) Dynamically ordered phase at 0=0.27. 



This symmetry breaking between the symmetric and asymmetric limit cycles has been the subject of intensive research 
over the last decade. It was first observed during numerical integration of a mean-field equation of motion for 
the magnetization of a ferromagnet in an oscillating field [^ , p2| . Since then, it has been observed and studied in 
numerous Monte Carlo (MC) simulations of kinetic Ising systems [ p^ , p7| -|33[| , as well as in further mean-field studies 
[p3| , p9 31 3^,0. It may also have been experimentally observed in ultrathin films of Co on Cu(OOl) |l2|. The 



results of these studies suggest that this symmetry 
phase transition. For recent reviews see Refs. |35|,|36f| 



breaking corresponds to a genuine continuous non-equilibrium 
Associated with the transition is a divergent time scale (critical 



slowing down) [^ and, for spatially extended systems, a divergent correlation length [ p7|j28| . Estimates for the critical 
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exponents and the universality class of the transition have recently become available | |27| , |28|,|37i after the successful 
application of finite-size scaling techniques borrowed from equilibrium critical phenomena ||38|-pl[] . 

The purpose of the present paper is to extend preliminary results and to provide more accurate estimates of 
the exponents for two-dimensional kinetic Ising systems in a square-wave oscillating field. The use of the square- 
wave field tests the universality of the dynamic phase transition (DPT) |2^,|2^, and it also significantly increases 
computational speed, compared to the more commonly used sinusoidal field. We further explore the universal aspects 
of the transition by varying the temperature and field amplitude within the multi-droplet regime, and we study the 
cross-over to the strong-field regime where the transition disappears. In obtaining our results, we rely on dynamic 
MC simulations. Computational methods are always helpful, especially when theoretical ideas are largely missing. 
There are cases, however, when even the use of standard equilibrium techniques, such as finite-size scaling requires 
some insight and building analogies between equilibrium and non-equilibrium systems [ p7| , p8[ . This is the case for 
our present study. No effective "Hamiltonian" was known before the completion of this work for the dynamic order- 
parameter (in the coarse-grained sense), from which the long-distance behavior of the model could be derived. This 
is a typical difficulty when dealing with systems far from equilibrium [^,^. Recently, however, a coarse-grained 
Hamiltonian has been derived for the dynamic order-parameter, supporting our results for the DPT. Similar to 
the previous work for sinusoidly oscillating fields p7| , p8t , we perform large-scale simulations and finite-size scaling to 
investigate the universal properties of the DPT. 

The remainder of the paper is organized as follows. In Sec. II we define the model and the observables of interest. 
Section III contains the Monte Carlo results and analyses. Conclusions and outlook are given in Sec. IV. 



II. MODEL AND RELEVANT OBSERVABLES 

To model spatially extended bistable systems in two dimensions, we study a nearest-neighbor kinetic Ising ferro- 
magnet on a Lxi square lattice with periodic boundary conditions. The model is defined by the Hamiltonian 

n^-JY,^^Sl-H{t)Y,S^, (1) 

where Si—±1 is the state of the zth spin, J > is the ferromagnetic interaction, runs over all nearest-neighbor 

pairs, runs over all lattice sites, and H(t) is an oscillating, spatially uniform applied field. The magnetization 
per site, 

"^w = iiE^'W' (2) 

2 = 1 

is the density conjugate to H{t). The temperature T is fixed below its zero-field critical value Tc (J/A:B7c=ln(l + 
V2)/2 where fee is Boltzmann's constant), so that the magnetization for H—0 has two degenerate spontaneous 
equilibrium values, ±r7isp(r). For nonzero fields the equilibrium magnetization has the same sign as while for H 
not too strong, the opposite magnetization direction is metastable and decays slowly towards equilibrium with time, 
as described in Sec. I. 

The dynamic used in this study, as well as in Refs. ||2^,^, is the Glauber single-spin-flip MC algorithm with updates 
at randomly chosen sites. Note that the random sequential update scheme corresponds to independent Poisson arrivals 
for the update attempts (discrete events) at each site. Thus, the arrival pattern is strongly asynchronous. The time 
unit is one MC step per spin (MCSS). Each attempted spin flip from Si to —Si is accepted with probability 

w( \ exp(-/3Ag,) 

Wis, -Si) = — . 3 

1 + exp(— pAi^i) 

Here AEi is the energy change resulting from acceptance, and /3= l/k-oT. For the largest system studied (L=512) 
we used a scalable massively parallel implementation of the algorithm for this asynchronous dynamics pq-Efl]. The 
parallel discrete-event scheme ensures that the underlying dynamic is not changed (that is, the update attempts 
are identical, independent Poisson arrivals at each site), while a substantial amount of parallelism is exploited. The 
parallel implementation was carried out on a Cray T3E, employing up to 256 processing elements. 
The dynamic order parameter is the period-averaged magnetization Ell], 
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where ti/2 is the half-period of the oscillating field, and the beginning of the period is chosen at a time when H{t) 
changes sign. Although the phase of the field does not influence the results reported in this paper, the choice made 
here is convenient in studies of the hysteresis loop-area distributions and consistent with Refs. [p5|-p^]. Analogously 
we also define the local order parameter 



1 



2ii/2 



S^it)dt , (5) 



which is the period-averaged spin at site i. For slowly varying fields the probability distribution of Q is sharply peaked 
at zero p7| , p8|| . We shall refer to this as the dynamically disordered phase. It is illustrated by the evolution of the 
magnetization in Fig. ||(a) and by the QwO time series in Fig. |^. For rapidly oscillating fields the distribution of Q 
becomes bimodal with two sharp peaks near ±msp(r), corresponding to the broken symmetry of the hysteresis loops 
p7| , p8[ . We shall refer to this as the dynamically ordered phase. It is illustrated in Fig. |^(b) and by the Q~TOsp~0(l) 
time series in Fig. ^ Near the DPT we use finite-size scaling analysis of MC data to estimate the critical exponents 
that characterize the transition. We also keep track of the normalized period-averaged internal energy (in units of J) 

13' 

i?=--^ i'-^^s,(t)s,(i)dt, (6) 

since it also exhibits important characteristics of the DPT. 

Previous studies of the DPT have used an applied field which varies sinusoidally in time. While sinusoidal or 
linear saw-tooth fields are the most common in experiments and are necessary to obtain a vanishing loop area in the 
low- frequency limit ||2^-|27[|, the wave form of the field should not affect universal aspects of the DPT. This should be 
so because the transition essentially depends on the competition between two time scales: the half-period ti/2 of the 
applied field, and the average time it takes the system to leave the metastable region near one of its two degenerate 
zero-field equilibria when a field of magnitude Hq and sign opposite to the magnetization is applied. This metastable 
lifetime, {t{T, Hq)), is estimated as the average first-passage time to zero magnetization. In the present paper we use 
a square-wave field of amplitude Hq. This has significant computational advantages over the sinusoidal field variation 
since we can use two look-up tables to determine the acceptance probabilities: one for H = +Hq and one for H = —Ho. 

In terms of the dimensionless half-period, 

e^t,/.2/{T{T,Ho)) , (7) 

the DPT should occur at a critical value 8c of order unity. Although Q can be changed by varying either ti/2, Hq, 
or T, in a first approximation we expect Oc to depend only weakly on Hq and T. This expectation will be confirmed 
in Sec. IV by simulations carried out at several values of Hq and T for different system sizes. 

In many studies of the DPT the transition has been approached by changing Hq or T [ p^ , p9|j3^j3^ . While this is 
correct in principle, (t(T, Hq)) depends strongly and nonlinearly on its arguments [p^ . We therefore prefer changing 
ti/2 at constant Hq and T |27 2^, as this gives more precise control over the distance from the transition. 



We focus on systems which are not only larger than the critical droplet, but also significantly larger than the typical 
droplet separation |15|. In this regime many supercritical droplets form and contribute to the decay of the metastable 
phase (the KJMA or Avrami theory for homogeneous nucleation [p"7|jl^ ]), as seen in Fig. |^. This is the only regime 
where the DPT is expected to exist. For small systems one observes subtle finite-size effects, not related to the DPT 
but rather to the stochastic single-droplet decay mod e |r^ ,|50|]. In the single-droplet regime, subject to a periodic 
applied field, the system exhibits stochastic resonance 



III. SIMULATION RESULTS 



A. Signs of the dynamic phase transition 



We performed extensive simulations and finite-size scaling analysis of the data on square lattices with L between 64 
and 512 at T=0.8Tc and Hq=0.3J. We also investigated the universality of the DPT within the multi-droplet regime 
for various fields and temperatures below the equilibrium critical temperature, using smaller systems with L from 64 
to 128. Typical runs near the DPT consist of 2x10^ full periods. For example, at T—0.8Tc and Hq—0.3J, where the 
critical half-period of the field is about 70 MCSS, this corresponds to 2.8x10'' MCSS. Away from the transition point, 
an order of magnitude shorter runs were sufficient to obtain high-quality statistics. 
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The system was initialized with all spins up and the square-wave external field started with the half-period in 
which H=—Hq. After some relaxation the system magnetization would reach a limit cycle [Fig. |^ (except for thermal 
fluctuations). In other words, Q [Eq. (^] (together with other period-averaged quantities) becomes a stationary 
stochastic process [Fig. |^. We discarded the first 1000 periods of the time series to exclude transients from the 
stationary-state averages. 

For large half-periods (0 ^ 6c) the magnetization switches every half-period [Fig. ^(a)] and QfnO, while for small 
half-periods (0 ^ &c) the magnetization does not have time to switch during a single half-period [Fig. ||(b)], resulting 
in |(5|~msp, as can be seen from the time series in Fig. The transition between the high- and low-frequency regimes 
is characterized by large fiuctuations in Q near 8c [Fig. |^ . 




_1.0 I • • • • • 1 

500 1000 1500 2000 
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FIG. 3. Time series of the order parameter Q &t T = O.STc and Ho — 0.3J for L = 128. Horizontal trace near Q = -1-1: 
= 0.27 < 0c [dynamically ordered phase, corresponding to Fig. ^b)]. Strongly fluctuating trace: = 0.98 ~ 0c (near the 
DPT). Horizontal trace near Q = 0: = 2.7 > 0c [dynamically disordered phase, corresponding to Fig. |^(a)]. 




FIG. 4. Configurations of the local order parameter {Q^} at T=0.8Tc and Ho=0.3J for L=128. (a) = 0.27 < 0c (dynami- 
cally ordered phase), (b) = O.98«0c (near the DPT), (c) = 2.7 > 0c (dynamically disordered phase). On the gray-scale 
black (white) corresponds to —1 (+1). 

To illustrate the spatial aspects of the transition we also show configurations of the local order parameter {Qi} in 
Fig. ^. Below 8c [Fig. ^(a)] the majority of spins spend most of their time in the state, i.e., in the metastable phase 
during the first half-period, and in the stable equilibrium phase during the second half-period (except for equilibrium 
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fluctuations). Thus, most of the Qi^+1. Droplets of Si=—1 that nucleate during the negative half-period and then 
decay back to +1 during the positive half-period show up as roughly circular gray spots in the figure. Since the 
spins near the center of such a droplet become negative first and revert to positive last, these spots appear darkest 
in the middle. Also, for not too large lattices one occasionally observes the full reversal of an ordered configuration 
{Qi}— *■{— Qi}, typical of finite, spatially extended systems undergoing symmetry breaking. Above Oc [Fig- H(c)] the 
system follows the field in every half-period (with some phase lag) and Qi~0 at all sites i. Near Oc [Fig. ^b)] there 
are large clusters of both Qi^-f 1 and 1 separated by "interfaces" where Qi«0. These large-scale structures 

remain reasonably stationary over several periods. 

For finite systems in the dynamically ordered phase the probability density of Q becomes bimodal. Thus, to capture 
symmetry breaking, one has to measure the average norm of Q as the order parameter, i.e., |3^. Figure ||(a) 

shows that this order parameter is of order unity for Q < 0^ and vanishes for O > Q^, except for finite-size effects. 
To characterize and quantify this transition in terms of critical exponents we employ the well-known technique of 
finite-size scaling [^-^. The quantity analogous to the susceptibility is the scaled variance of the dynamic order 
parameter p7|j2^ , 

= L' {{Q')l - ml) . (8) 

Note that for our system the field conjugate to Q and a corresponding fiuctuation-dissipation theorem are not known, 
hence we cannot measure the susceptibility directly. For finite systems Xl has a characteristic peak near Qc [see 
Fig. ^(b)] which increases in height with increasing L, while no finite-size effects can be observed for Q ^ and 
O ^ Oc. This implies the existence of a divergent length scale, possibly the correlation length which governs the 
long-distance behavior of the local order-parameter correlations {QiQj). The location of the maximum in X'^ also 
shifts with L, which gives further important information about the critical exponents. 
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FIG. 5. Finite-size behavior of the order parameter at T^O.STc and //o=0.3J for various system sizes, (a) The order 
parameter {|Q|)l. (b) The scaled variance of the order parameter, X^, as defined in Eq. (^. (c) Same as (b) on lin-log scale 
to provide an enhanced view of the peaks for smaller systems. 



The normalized stationary time-displaced autocorrelation function of the order parameter, 



.Q, ^ ^ {Q{m{j+n))-{Q{3)? 



(9) 



provides further insights into the DPT as the system exhibits critical slowing down [Fig. ^. This can be seen as 
increasing correlation times with increasing system sizes. In Sec. III.D we provide a quantitative analysis of the 
correlation times. 

We also measured the period-averaged internal energy [Eq. (||)] and its fluctuations [^ 



XE = L^{E^),^{E)l) 



(10) 



as can be seen in Fig. 0. The peaks of these fluctuations exhibit a slow increase with the system size (compared to 
the order-parameter fluctuations) , as one may anticipate by analogy with the equilibrium heat capacity. 
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FIG. 6. Critical slowing down for the order parameter at T=0.8rc and i/o=0.3J at 0=Oc, as shown by the normalized 
autocorrelation function C^(n). 
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FIG. 7. Finite-size behavior of the period-averaged internal energy at T^O.STc and Ho=0.3J for various system sizes, (a) 
The period-averaged internal energy, (b) The scaled energy variance, Xj^ as defined in Eq. (^o|). 



B. Finite-size scaling 



Scaling laws and finite-size scaling for equilibrium systems with an a-priori known Hamiltonian can be systematically 
derived using the concepts of the free energy and the renormalization group The kinetic Ising model with 

the explicitly time-dependent Hamiltonian, Eq. (|^), is driven far from equilibrium. Although the order-parameter 
distribution P{Q) is stationary, the effective Hamiltonian controlling its fixed-point behavior has not been known 
until recently (after the completion of this study). Motivated by the similarity of the finite-size effects shown 
in Figs. to those characteristic of a typical continuous phase transition, we borrow the corresponding scaling 
assumptions from equilibrium finite-size scaling. For our model the quantity analogous to the reduced temperature 
in equilibrium systems (i.e., the distance from the infinite-system critical point) is 

„ ie-e,i , , 

(11) 



Finite-size scaling theory provides simple scaling relations for the observables for finite systems in the critical regime 

{\Q\)L=L-f'/-T±{eL^/'') (12) 
Xl ^ L^'^g^iOL^'") (13) 
Xl = c,\n[LJ±{6L^''')) , (14) 
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where T±^ Q±, and J7± are scaling functions and the + (— ) index refers to > 0c (0 < 0c)- The logarithmic 
scaling in is motivated by the very slow divergence of the scaled period-averaged energy variance [Fig. 0(b)]. The 
above formulation of scaling is explicitly based on the infinite-system critical point 0c , which can be estimated with 
far greater accuracy than the location of the maximum of the order-parameter fluctuations for the individual finite 
system sizes. We use the fourth-order cuniulant intersection method [ ^J40| to estimate the value of ©c at which the 
transition occurs in an infinite system. In order to do this, we plot 

as a function of for several system sizes as shown in Fig. ^. For the largest system {L—512) the statistical uncertainty 
in Ul was too large to use it to obtain estimates for the crossing. Our estimate for the dimcnsionlcss critical half- 
period, based on the remaining five system sizes, is 0c=O.918±O.OO5 with a fixed-point value ?7*=0.611±0.003 for the 
cumulant [Fig^(b)]. 
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FIG. 8. (a) The fourth-order cumulant as defined in Eq. ( p^ at T^O.STc and J/o=0.3J for various system sizes, (b) The 
region around the cumulant crossing in (a) enlarged. The horizontal and vertical dashed lines indicate the fixed-point value 
i7*=0.611 and the scaled critical half-period Oc=0.918, respectively. 



Then at 0c the scaling forms Eqs. ( p^[l^ yield 



Xl (XC2 + ci ln(L) , 



(16) 
(17) 
(18) 



which enable us to estimate the exponent ratios P/v and 7/1^, and to directly check the postulated logarithmic 
divergence in the period-averaged energy fluctuations. Plotting {\Q\)l and X'^ at ©c and utilizing a weighted hnear 
least-squares fit to the logarithmic data yields /3/j/=0.126 ± 0.005 [Fig. |9|(a)] and 7/1^=1. 74 ± 0.05 [Fig. ||(b)]. Note 
that these values are extremely close (within statistical errors) to the corresponding ratios for the equilibrium two- 
dimensional Ising universality class, /?/i/=l/8=0.125 and 7/j/=7/4=1.75. Further, the straight line in Fig. ^(c) 
indicates the slow logarithmic divergence of Xf^ at the critical point. In addition to the scaling at ©c, we also checked 
the divergences of the peaks of the fluctuations, (X^)poak and (Xf')pcak, since they asymptotically should follow the 
same scaling laws, Eqs. ( p7| ) and (p^), respectively. The measured exponent 7/1^=1. 78 ± 0.05 for (^^)peak and the 
logarithmic divergence for (Xf')poak agree to within the statistical errors with the results obtained at ©c, as can be 
seen in Fig. ^b) and (c), respectively. 

Prom the finite-system shifting of the transition one can estimate the correlation- length exponent u by tracking the 
shift in the location of the maximum in X'^: 



|©c(L) - ©c| oc 



-l/v 



(19) 



where ©c(^) is the location of the peak for finite systems. However, the precision of this method for our data is very 
poor, due to limited resolution in finding the locations of the maxima and consequently the large relative errors in 
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|8c(i) — 6c|- Excluding the smallest (due to strong corrections to scaling) and the largest systems (due to very poor 
resolution and extremely large statistical error), we obtain i/=0.87±0.4, but the large error estimate obviously implies 
rather poor accuracy [Fig. |0[. 
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FIG. 9. Critical exponent estimates at T—0.8Tc and Hq—O.SJ. Straight lines are the weighted least-square fits, (a) De- 
termining f3/u through the finite-size effects of the order parameter, based on Eq. ( p^ (log-log plot), (b) Determining j/iy 
through the finite-size effects of the order-parameter fiuctuations, based on Eq. (|l^) (log-log plot), (c) Showing the logarithmic 
divergence of the period-averaged energy fluctuations, based on Eq. (jl^) (log-lin plot). 
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FIG. 10. Exponent estimate for u at T—0.8Tc and 1/0=0.3 J, based on Eq. ( pj| ) (log- log plot). The dot-dashed line is a 
weighted least-squares fit (excluding the smallest and the largest system) yielding the slope l/!/=1.15 ± 0.6 (^'=0.87 ± 0.4). 
The dashed line represents the "optimal" value for this exponent, using the best quality data collapse for the scaling function 
Eq. ( p^ as discussed in the text, yielding !/«/=!. 05 (i^=0.95). The solid line represents the two-dimensional equilibrium Ising 
exponent i^=1.0. 

To obtain a more complete picture of how well the scaling relations in Eqs. jl^ ) and ( p^ hold, we plot {\Q\)lL'^^'^ 
[Fig. |l|(a)] and X^L'''/'' [Fig. |ll|(b)] vs 91^/" [1^. For the exponent ratios we used /3/j/=l/8=0.125, and 
7/^=7/4=1.75, since our estimate for those (within small statistical errors) implied that they take on the equi- 
librium two-dimensional Ising universal values. Most importantly, we used various values of v between 0.5 and 1.2 
to find the best data collapse as observed visually, since our estimate for this exponent was far from reliable. The 
"optimal" value obtained this way (by showing scaling plots to group members who did not know the particular values 
of v used), and used in Fig. |ll|(a) and (b), is 1^=0.95 ± 0.15. Full scaling plots using the exact Ising exponents are 
also shown in Fig. |l^(c) and (d), and they result in similarly good data collapse. 
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FIG. 11. Finite-size scaling (full data collapse) at T—0.8Tc and Ho~0.3J using f3/v=l/S, y/f—J/'l (two-dimensional equi- 
librium Ising values), and v—O.Qb (which yields the best quality data collapse), (a) For the order parameter (|Q|)l (log-log 
plot), (b) For the scaled order-parameter variance (log-log plot). Straight lines in both graphs represent the asymptotic 
large- argument behaviors of the scaling functions J-± and Q± given in Eqs. ( |l2[ ) and (^3|), respectively. Figures (c) and (d) are 
the same as (a) and (b), except that the exact Ising exponent v=1.0 is used. 



C. Order-parameter histograms at criticality 

We devote this subsection to analyzing the universal characteristics of the full order-parameter distribution, P{Q), 
at the critical point. This distribution is bimodal for finite systems if observed for sufficiently long times (Fig. Il2) 
It is more convenient to focus on the distribution of \Q\, avoiding the effect of the insufficient number of 
switching events between the two symmetry-broken phases for large systems, which causes the skewness in Fig. [l^(b). 
Figure E3(a) shows the order-parameter distributions -Pl (|Q|) a-t the critical point for various system sizes. Finite-size 
scaling arguments suggest that at 9c 

PlUQD ^ L^/'^ViLf'/'^m . (20) 

Thus, the scaled distributions, L-'^^'' Pl{\Q\) vs x=\Q\L'^^'', should faU on the same curve V(x) for different system 
sizes. Again, we used P/v = 1/8. The quality of the data collapse is quite impressive [Fig. |l^(b)], with deviations 
mainly observed for the smallest L and the largest values of \Q\ [Fig. |l^(c)], possibly as a result of corrections to 
scaling. 
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FIG. 12. (a) Short segment of the order-parameter time series at T=0.8rc and Ho=0-SJ for an L=128 system at | p2[ . 
(b) Order-parameter histogram for the same parameters. 

What we find somewhat surprising, is that the distribution appears to be identical (except for stronger corrections to 
scahng at the DPT) to that of the equihbrium two-dimensional Ising model on a square lattice with periodic boundary 
conditions at criticality, without a need for any additional scaling parameters. We checked this by performing standard 
equilibrium two-dimensional Ising simulations with Glauber dynamics and system sizes ranging from L=64 to L—128, 
and also by comparing our scaled DPT order-parameter histograms to the high-precision two-dimensional equilibrium 
Ising MC data of Ref. |Q (Fig- ^H)- We had expected the shapes of the distributions to be identical for the DPT and 
equilibrium Ising model, as a consequence of the identical values for the cumulant fixed-point value U*. However, it 
is not obvious to us why the microscopic length scales in the DPT and the equilibrium Ising model also appear to be 
identical, as evidenced by the absence of the need for an additional scaling, L^L/a (a being the microscopic length 
scale in the DPT). 




FIG. 13. Order-parameter histograms Pl(\Q\) at T^O.STc and Ho=Q.?>J for various system sizes at Oc. The thin solid lines 
represent two-dimensional equilibnum Ising order-parameter histograms at the critical point for L—64, 90, and 128 on all three 
graphs, (a) Order-parameter distributions, (b) Scaled order-parameter distributions, according to Eq. (|^ with /3/i/=l/8. 
The bold solid line is the corresponding (Monte Carlo) two-dimensional equilibrium Ising distribution without any additional 
scaling parameters |5^. (c) Same as (b) on lin-log scales to enhance the view of the corrections to scaling for small systems at 
large \Q\. 



D. Critical slowing down 

Computing the stationary autocorrelation function given by Eq. (^) , we already pointed out that at Qc the correla- 
tion time increases fast with system size [Fig. ^. Correlation times are typically extracted from an exponential decay 
as 

C^(n)(xe-"/--^ (21) 
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and they are expected to be finite for finite systems. The correlation time is also well defined in the L— >oo limit 
away from the transition. However, it diverges with L at the transition point as 

(X L% (22) 

where z is the dynamical critical exponent. For not too late times we had reasonable statistics including the larger 
systems (up to L—256) to fit the usual exponential decay [Fig. |l^(a)]. Then plotting the correlation times vs L 
yields the dynamic exponent z = 1.91 ± 0.15, as shown in Fig. |14|(b). This value is within two standard deviations of 
most estimates for the dynamic exponent of the two-dimensional equilibrium Ising model with local dynamics [p4[ . 




80 100 120 ,140 

n [periods 

FIG. 14. Critical slowing down for the order parameter at T—0.8Tc and Hq—OSJ at Q=Qc. (a) The normahzed autocor- 
relation function on lin-log scale for early times. The straight lines are fits to exponential decays according to Eq. (pl|). (b) 
Determining the dynamic exponent, z, using a power-law fit to Eq. ( [22[ ) (log-log plot). 



E. Universality for various temperatures and fields and cross-over to tlie strong-field regime 

The underlying ingredient for the spatially extended bistable systems exhibiting a DPT is the local metastability 
(and the corresponding characteristic time spent in the metastable "free-energy well" ) in the presence of an external 
field. This, in turn, provides a competition between time scales if the system is driven by a periodic field. Based 
on this, we expect that sufficiently large systems (in which many droplets contribute to the decay of the metastable 
phase) exhibit the DPT at a half-period ii/2 comparable to the metastable lifetime {t{T, Hq)). In other words, we 
expect the critical dimensionless half-period to be of order one, 8c ~ 0{l). 

To test this expectation, we performed simulations at T—0.8Tc for field amplitudes ranging from 0.3J to infinity 
with system sizes L=64, 90, and 128. [Ho=oo corresponds to the Glauber spin-flip probabilities Eq. (^) being equal 
to (1) depending on whether the spin is parallel (anti parallel) to the external field, with no influence from the 
configuration of the neighboring spins.] We further performed runs at r=0.9Tc, T=0.6Tc, and r=0.5rc for various 
field amplitudes and system sizes i = 64 and 90. The typical run length was 2x10^ periods. The purpose of these 
runs was to explore the universal nature of the DPT in the multi-droplet regime, and the crossover to the strong-field 
regime where the DPT should disappear. 

In the strong-field regime the droplet picture breaks down since the individual spins are decoupled. Thus, the 
metastable phase no longer exists, and the decay of the phase having opposite sign to the external field approaches 
a simple exponential form (which becomes exact in the Hq—^oo limit). In the Appendix we show that under these 
conditions the system magnetization always relaxes to a symmetric limit cycle with Q^O for all frequencies, thus, no 
DPT can exist. 

Figures |l^(a), (b), and (c) show the order parameter vs the dimensionless half-period for L — 6A and a range of 
field amplitudes at T=0.8Tc, T=0.6Tc, and T—Q.5Tc, respectively. The typical order-parameter profile where the 
system exhibits the DPT prevails up to some temperature dependent cross-over field amplitude (T) [filled symbols 
in Fig. |l^. For Hq > Hx (T) the underlying decay mechanism belongs to the strong-field regime, and correspondingly 
the DPT disappears as expected. 
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FIG. 15. DPT in the multi-droplet regime and strong-field cross-over for various temperatures and field amplitudes with 
I/=64. Order-parameter profiles with filled symbols exhibit a DPT. Corresponding lifetimes in units of MCSS are also indicated 
in the figures, (a) T^O.STc (b) T=0.&Tc. (c) r=0.5rc. 



One can trace this crossover from the multi-droplet to the strong-field regime by plotting 0c vs Ht^/ J [Fig. p^(a)] 
and vs (r) [Fig. |l^(b)] for fixed temperatures. Here again we employed the fourth-order cumulant intersection method 
using L—64 and L=90 to identify the infinite-system transition point, 0c. The crossover to the strong-field regime is 
indicated by the drop in 0c for large fields (small lifetimes). 




H/J <%> T/T^. 



FIG. 16. Scaled critical half-period in the multi-droplet /DPT regime: (a) as a function of the external field amplitude Hq, 
(b) as a function of the lifetime (r). In (a) the corresponding empty symbols for each temperature on the horizontal axis 
represent the points which are the MC upper bounds for the cross-over field amplitude, Hx{T), beyond which no DPT can 
be found for any non-zero 0. In (b) the corresponding empty symbols are the lower bound for the lifetime below which no 
DPT is observed, (c) Metastable phase diagram for the kinetic Ising model. The stars connected by a dashed line give a 
numerical estimate for the dynamic spinodal for our smallest system, L=64, that separates the multi-droplet (MD) from the 
single-droplet (SD) regime |13]. The solid line is an analytic estimate for the "mean-field spinodal" which separates the MD 
from the strong-field (SF) regime The circles represent the temperature and field amplitude values at which we ran our 

simulations. The filled circles indicate points where the system exhibits a DPT, and empty circles represent points where it 
does not. 



The DPT in the multidroplet regime {Ho<Hx{T)) appears to be universal, as can be seen from the scaled order- 



parameter plot in Fig. 17, Note that both graphs contain 28 DPT data sets: three different temperatures, with four 
different field amplitudes for each, and at least two different system sizes (L = 64, 90, and 128 at r=0.8Tc, and 
L = 64 and 90 at r=0.6Tc and T=0.5Tc) for all these parameters! While the slopes in the asymptotic scaling regime 
appear to be the same, the small parallel shift may be the result of the non-universal critical amplitudes at different 
temperatures and fields, or simply our inaccuracy in determining Oc(T, Hq), due to the relatively short runs. 
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FIG. 17. Universality of the DPT in the multi-droplet regime, (a) Finite-size scaling (full data collapse) for the order 
parameter (|(5|)l (log-log plot). The figure contains 28 DPT data sets: three different temperatures, with four different field 
amplitudes for each, and at least two different system sizes (L = 64, 90, and 128 at r=0.8rc, and L = 64 and 90 at r=0.6rc and 
T=0.5Tc) for all these parameters. The temperature and field values were chosen such that the system is in the multi-droplet 
(MD) regime [filled circles in Fig. ^6|(c)]. We used /3/!/=l/8, 7/i/=7/4 (two-dimensional equilibrium Ising values), and !/=0.95, 
the optimal value, as described in the Section III.B. Straight lines represent the asymptotic large- argument behaviors of the 
scaling functions T± given by Eq. (|l2). Figure (b) is the same as (a), except that the exact Ising exponent v=l.Q is used. 



On the other hand, as expected, the system shows no singularity and (IQI)^— — >°0 for any non-zero frequency in 
the strong-field regime, Ho>Hx{T) (see the Appendix). Here, the finite-size eff'ects simply reflect the central-limit 
theorem, i.e., {Q^)^0{l/L^), or (|Q|)~0(1/L). Figures [l8|(a) and (b) illustrate this at T=0.8Tc, Ho^l.bJ and at 
T=0.8Tc, Hq=oo, respectively. One also expects that a similar behavior prevails at T>Tc for any field amplitude, i.e., 
i?x(r) vanishes for T>Tc. We checked this for T=l.irc and iJo=0.05J, and the results confirm the {\Q\)^0{l/L) 
scaling of the order parameter for all frequencies, showing none of the characteristic finite-size effects of a DPT 
[Fig. |l^(c)]. This result is in distinct disagreement with older work such as Rcfs. |^,^, which report observations 
of the DPT at temperatures considerably above Tc- For smaller system sizes, however, one may observe nontrivial 
resonance Ipal. 
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FIG. 18. Order parameter (|(5|)l behavior outside the multi-droplet regime, (a) Strong-field behavior at r=0.8rc, Ho=1.5J 
({r)=4.3 MCSS). (b) Strong-field behavior at r=0.8rc, i^o^oo ((r)=ln2 MCSS). (c) High-temperature behavior at T=l.lTc, 
-^0=0. 05 ({r)=81.8 MCSS). The insets in all three graphs show data collapse for the scaled order parameter {\Q\)lL. 



Due to the expected complications of a divergent equilibrium correlation length, we did not perform simulations at 
Tc. However, we conjecture that Hx{T) vanishes at Tc- We expect this conjecture to be extremely difficult to prove 
or disprove numerically, due to very large and possibly complicated finite-size effects. 
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IV. CONCLUSIONS AND OUTLOOK 



In this paper we have studied the hysteretic response of a spatially extended bistable system exhibiting a DPT. 
Our model system is the two-dimensional kinetic Ising ferromagnet below its equilibrium critical temperature subject 
to a periodic square-wave applied field. The results indicate that for field amplitudes and temperatures such that 
the metastable phase decays via the multi-droplet mechanism, the system undergoes a continuous dynamic phase 
transition when the half-period of the field, ti/2, is comparable to the metastable lifetime, {t{T, Hq)). Thus, the 
critical value 6c of the dimensionless half-period defined in Eq. is of order unity. As 6 is increased beyond 9c, the 
order parameter {\Q\) (the expectation value of the norm of the period-averaged magnetization) vanishes [Fig. ||(a)], 
displaying singular behavior at the critical point, as shown in Fig. |^(b) and (c). 

The characteristic finite-size effects in the order parameter and its fluctuations indicate that there is a divergent 
correlation length associated with the transition. We used standard finite-size scaling techniques adopted from the 
theory of equilibrium phase transitions. We estimated Oc and the critical exponents /3, 7, and v from relatively 
high precision data for system sizes between L = 64 and 512 at T=0.8Tc and Ho=0.3J. Our best estimates are 
(3 lv=Q.\2Q ± 0.005, 7/1^=1. 74 ± 0.05, and ^=0.95 ± 0.15. These values agree within statistical errors with those 
previously obtained with a sinusoidally oscillating field ||2^,|2^ , providing strong evidence that the shape of the field 
oscillation does not affect the universal aspects of the DPT. Observing the stationary autocorrelation function we 
also saw that at the transition point the system exhibits critical slowing down governed by the dynamic exponent 
2;=1.91±0.15. This is also very close to the corresponding exponent z=2.12(5), measured in standard two-dimensional 
Ising simulations with local dynamics jEJ]. Our best values for the exponent ratios P/v and ^/v are given with 
relatively high confidence, while for u it is rather poor. In this sense tracking down the exponent v and obtaining an 
accurate estimate for it remains elusive. Note, however, that we could only rely on the standard (single spin-flip) MC 
algorithm, since we had to preserve the underlying dynamics. Using more sophisticated algorithms to avoid critical 
slowing down as seen in Figs. I and |l|, would require an underlying "Hamiltonian" for the corresponding local order 
parameter {Qi}, which is not yet known. While in a coarse-grained/universal sense a 0^ Hamiltonian is supported 
by our data, it does not point to any one particular microscopic Hamiltonian for the microscopic order parameter 
{Qi}. However, a 0^ coarse-grained Hamiltonian for \Qi \ has been recently derived starting from the time-dependent 
Ginzburg-Landau equation for the magnetization 

Of the known universality classes, our exponent estimates for the DPT are closest (and within the statistical errors) 
to those of the the two-dimensional equilibrium Ising model: /3/i^=l/8=0.125, 7/z^=7/4=1.75, and v=\. Consequently, 
our measured exponent ratios satisfy the hyperscaling relation 

2(/3/i/) -1-7/1/ = 1.99 ±0.05 w d , (23) 

where d=2 is the spatial dimension [Q. Further, the fixed-point value of the fourth-order cumulant, [/*=0.611±0.003, 
is also extremely close to that of the Ising model, t/* =0.6106901 (5) Q. These findings provide conclusive evidence 
that the DPT indeed corresponds to a non-trivial fixed point. We tested the full data collapse for the scaled order 
parameter and its variance, as shown in Fig. and it confirmed the existence of the universal scaling functions 
given by Eqs. (|l2) and (|l^). Also, at the critical frequency, the order-parameter distributions follow finite-size scaling 



predictions, Eq. (2^), as shown in Fig. |T^. More surprisingly, the critical DPT order-parameter distributions fall on 
that of the two-dimensional equilibrium Ising model at the critical temperature (except for stronger corrections to 
scaling), without any additional fitting of the underlying microscopic length scale. 

While our finite-size scaling data clearly indicate the existence of a divergent length scale, we did not measure the 
correlation length for the local order parameter directly. Future studies may include extracting the correlation length, 
fg, from the (QiQj) correlations (or from the corresponding structure factor). This approach would also provide 
another way to measure the exponent v by plotting vs 9 in the critical regime for large systems, and assuming 

We also studied the universal aspects of the DPT at other temperatures and fields. Shorter runs at r=0.8Tc, 
T=0.6Tc, and r=0.5Tc for field amplitudes Ho<Hy (T) also confirmed scaling and the universality of the DPT 
[Fig. |l^. The condition for the field amplitude implies that the system only exhibits a DPT in the multi-droplet 
regime. For Ho>Hx{T) strong-field behavior governs the decay of the magnetization, and the DPT disappears, as 
indicated by Fig. [l6| and Fig. ^(a),(b). We also found that the high-temperature phase is qualitatively similar to the 
strong-field regime in that there is no sign of a DPT for T>Tc [Fig. [l^(c)]. 

One may ask how general the phenomenon of a DPT is in spatially extended bistable systems, subject to a periodic 
applied "field" which drives the system between its metastable and stable "wells." It is possible that having up-down 
symmetry of the period-averaged "magnetization" is sufficient for possessing a Hamiltonian at the coarse-grained 
level, even if the system is driven far away from equilibrium and is microscopically irreversible |]58|| . Future research 
can address this question by studying other systems (not necessarily ferromagnets) that exhibit hysteresis. 



15 



ACKNOWLEDGMENTS 



We thank S. W. Sides, S. J. Mitchell, and G. Brown for stimulating discussions. We would like to thank W. 
Janke for providing us with data from Ref. p3|| for comparison of the critical order-parameter distributions. We 
acknowledge support by the US Department of Energy through the former Supercomputer Computations Research 
Institute, by the Center for Materials Research and Technology at Florida State University, and by the US National 
Science Foundation through Grants No. DMR-9634873, DMR-9871455, and DMR-9981815. This research also used 
resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science 
of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098. 



APPENDIX: LOW-FREQUENCY AND STRONG-FIELD MEAN-FIELD APPROXIMATION 

For simplicity, in the following we assume that the magnetization decays from m—+l to m——l after a single field 
reversal {Hq—>-—Ho). This is a good approximation below the equilibrium critical temperature (rrispssl), and for any 
temperature when Hq—i-oo. Further, we assume that the volume fraction of mcta- or unstable spins follows a simple 
monotonic decay, ip{t). 

In terms of the volume fraction of positive spins, (j){t), the magnetization can be written as 

m{t) = 2(j){t) ~ 1 . (Al) 

Subject to a square-wave field, 

H{t) = h^° , (A2) 

in the first (second) half-period the volume fraction of the positive (negative) spins decays according to ip{t). Thus, 
in each period (measuring time t from the beginning of the period) 

' I 1 - [1 - 0(^1/2)] ^{t - <i/2) tl/2 < t < 2ii/2 • ^"^^^ 

Using this approximation, one directly obtains a linear mapping 

(l^n+l = 1 - [1 - (pn^{tl/2)Mh/2) , (A4) 

where (j)n=(j){2nti/2) is the volume fraction of the positive spins at the beginning of the nth period, n = 0, 1,2, . . .. 
The stationary value of this quantity is 

0* = lim (t)n = , A, r . (A5) 

Consequently, the magnetization reaches a stationary limit cycle 

m(<) « i ^ , . (A6) 

l+vk,,M t - *l/2) tl/2 < t < 2ti/2 

In this limit cycle the magnetization oscillates about zero, 

m(0) = -m{t,/,) = i-^^ (A7) 
1 + ^(^1/2) 

and the symmetry of the magnetization, m(t ± ii/2)=^"T'(^) implies 

Q = (p m{t)dt = . (A8) 

2^1/2 J 

This corresponds to the symmetric (dynamically disordered) phase. Note that this symmetric phase is always reached 
when (p decreases monotonically from unity at t=Q to zero as t—^oo 
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In the multi-droplet regime, the volume fraction of the metastable phase decays according to Avrami's law p7| , p^ . 
In the low-frequency limit, ti/23>(T) (03>1), each half-period almost always contains a complete metastable decay 
[Fig. 11(a)] . Avrami's law for the metastable volume fraction in each half-period can then be directly applied using 
(^(t)=(^^,(i)«e-(i"2)tV(^>' jl|Jl|]. This functional form for ip^sit) breaks down whe n ti/2 becomes comparable to 
(t) (0«1), thus this simple mean- field approximation cannot predict any instability related to the DPT. 

In the Hq^oo limit the individual spins become decoupled. Then one can obtain (^{t)=e^'-'"^-'*/^'^^ which is exact 
for all frequencies. Further, this exponential decay is a good approximation everywhere in the strong-field regime, 
thus, no DPT can exist there. 
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